Rt per provincia¶

Rt di COVID-19 stimato per le province italiane.

Max Pierini, Alessio Pamovio & NOTIZIÆ Telegram Channel

Method¶

A simple method is presented to estimate effective reproduction number $R_t$ of COVID-19 in italian provinces with a Markov chain Monte Carlo and Poisson likelihood parametrized on daily new cases.

Details¶

New cases $y_r$ for each $r$ italian province (source Dipartimento Protezione Civile), will be smoothed with rolling mean (gaussian, window 7, std 2). Smoothed new cases will be adjusted to be $>0$ to avoid negative values (due to data error corrections).

For each day $t$, smoothed new cases $y_{r,t}$ will be supposed distributed as Poisson with $\lambda_{r,t}$ parameter

$$ y_{r,t} \sim \mathcal{P}(\lambda_{r,t}) $$

where $\lambda_{r,t}$ is defined by the serial interval inverse $\gamma$, previous day smoothed new cases $k_{r,t-1}$ and effective reproduction number in time $R_{r,t}$ (ref: Bettencourt & Ribeiro 2008)

$$ \lambda_{r,t} = k_{r,t-1} e^{\gamma (R_{r,t} - 1)} $$

Parameters $R_{{p,t}}$, for each province $p$ and time $t$ (day) will be distributed as half normal with mean equal to previous day posteriors $R_{{p,t-1}}$ and precision $\tau$

$$ R_{{p,t}} \sim \mathcal{{N}}^+(R_{{p,t-1}}, \tau) $$

where, first day $R_{{p,0}}$ (outcome) is set to zero

$$ R_{{p,0}} = 0 $$

and $\tau$ is defined as $1 / \sigma^2$, where $\sigma$ is the standard deviation of $R_t$ priors, previuosly estimated for italian provinces (using the same method but with overarching uninformative prior distribution of $\tau$) as

$$ \sigma \sim \mathcal{N}(\mu_{=0.299}, \sigma_{=0.004}) $$

If previous new cases are zero $k_{p,t-1}=0$, parameter $R_{p,t}$ is undefined, given the chosen function for $\lambda_{p,t}$ parameter of Poisson likelihood, even if it should be $R_{p,t}=0$ (no new cases means null effective reproduction number). Thus, in these cases, priors of $R_{p,t}$ will be forced to

$$ R_{p,t} \sim \mathcal{N}^+(0, \tau) \;,\; k_{p,t-1}=0 $$

and previous new cases will be forced to $k_{p,t-1}=1$, so that $\lambda_{p,t}$ will be

$$ \lambda_{p,t} = e^{\gamma( \mathcal{N}^+(0, \tau) - 1 )} \;,\; k_{p,t-1}=0 $$

model {
    # Overarching Rt standard deviation
    sigma_R ~ dnorm( 0.29929894605552754 , 1 / 0.004218809765237365 )
    tau_R <- 1 / sigma_R^2


    # Serial interval distribution
    SI_mu <- 7.5
    SI_sd <- 3.4
    SI_sh <- SI_mu^2 / SI_sd^2
    SI_ra <- SI_mu / SI_sd^2
    SI ~ dgamma( SI_sh , SI_ra )
    gamma <- 1 / SI

    # First Rt prior
    R[1] <- 0
    for ( t in 2:T ) {
        # Rt prior for k>0
        Rpr[t] ~ dnorm( R[t-1] , tau_R )  T(0,)
        # Rt prior for k=0
        Rnu[t] ~ dnorm( 0 , tau_R )  T(0,)

        # Define Rt prior
        R[t] <- ifelse( k[t-1]==0 , Rnu[t] , Rpr[t] )
        # Avoid k=0 (undefined Rt)
        K[t] <- ifelse( k[t-1]==0, 1 , k[t-1] )

        # Poisson likelihood
        lambda[t] <- K[t] * exp( gamma * ( R[t] - 1 ) )
        y[t] ~ dpois( lambda[t] )
    }
}

model {
    # Overarching Rt standard deviation
    sigma_R ~ dnorm( 0.29929894605552754 , 1 / 0.004218809765237365 )
    tau_R <- 1 / sigma_R^2


    # Serial interval distribution
    SI_mu <- 7.5
    SI_sd <- 3.4
    SI_sh <- SI_mu^2 / SI_sd^2
    SI_ra <- SI_mu / SI_sd^2
    SI ~ dgamma( SI_sh , SI_ra )
    gamma <- 1 / SI

    # First Rt prior
    R[1] <- 0
    for ( t in 2:T ) {
        # Rt prior for k>0
        Rpr[t] ~ dnorm( R[t-1] , tau_R )  T(0,)
        # Rt prior for k=0
        Rnu[t] ~ dnorm( 0 , tau_R )  T(0,)

        # Define Rt prior
        R[t] <- ifelse( k[t-1]==0 , Rnu[t] , Rpr[t] )
        # Avoid k=0 (undefined Rt)
        K[t] <- ifelse( k[t-1]==0, 1 , k[t-1] )

        # Poisson likelihood
        lambda[t] <- K[t] * exp( gamma * ( R[t] - 1 ) )
        y[t] ~ dpois( lambda[t] )
    }
}

DONE!

Time series¶

Latest Rt (IQR)¶

Latest Rt (HPDI)¶

Latest Rt (MAP)¶

-> 2020-11-09 19:58:45.833174
<- 2020-11-09 21:13:03.605529
Completed in 1:14:17.772355